FlopyとPHREEQCで熱水地球化学シミュレーション #3:解析・可視化編

カップリングシミュレーションで得られたデータを読み解く。破過曲線・鉱物量の空間プロファイル・飽和指数の「意味」と「読み方」を、グラフの見方から丁寧に解説する。
Groundwater
Python
MODFLOW
Geochemistry
Simulation
地下水
作者

DeepFlow

公開

2026年4月5日

数字を「見える化」する意味

第2回でカップリングループを動かし、シミュレーションのデータが蓄積されました。しかし数値の羅列だけでは、カラムの中で何が起きているかは見えてきません。

地球化学シミュレーションの醍醐味は、反応フロントがカラムを移動する様子を時空間でとらえることです。

「反応フロント」とは、岩石と水の反応が最も活発に起きている境界面のこと。熱水がカラムに流れ込むにつれて、この境界が少しずつ奥へ移動していきます。それをグラフで「見える化」するのが今回のテーマです。

今回は4種類のグラフを取り上げます:

グラフの種類 横軸 縦軸 何がわかるか
破過曲線(pH) 時間 [h] pH 反応が出口に届くまでの「遅れ」
Ca²⁺ 濃度と飽和指数 時間 [h] Ca²⁺・SI 溶解がどのペースで進んでいるか
鉱物量の空間プロファイル カラム位置 [m] 方解石の残量 どこで溶解が最も進んでいるか
4パネルサマリー 時間・位置 各指標 全体像を一枚で俯瞰する

観測ポイントの設計:入口と出口を比べる

グラフを読む前に、どこを「観測ポイント」にするかを決めておきます。今回はカラムの両端、2点に注目します。

入口 [Cell 1]                        出口 [Cell 15]
    │                                     │
    ▼                                     ▼
新鮮な熱水が                        岩石と長時間
岩石に接触する                       反応した後の流体
「反応の最前線」                     「反応済みの流体」

この2点を重ねてグラフに描くことで、反応がカラムを通り抜けるまでの時間的な遅れ(ラグ) が一目でわかります。「入口ではもう反応が落ち着いているのに、出口にはまだ変化が届いていない」という状況が視覚的に確認できます。


グラフ 1:破過曲線(pH の時間変化)

破過曲線とは?

破過曲線(Breakthrough Curve) とは、出口での水質が時間とともにどう変化するかを示したグラフです。もともとは「注入した成分がカラムの出口に”破過”(突き抜けて)出てくる様子」を描いたものです。

横軸に時間、縦軸に pH をとり、入口(青の実線)と出口(アンバーの破線)を重ねて描きます。

読み方のポイント

① 定常状態への到達時刻

出口の pH が時間とともに変化し続けているうちは、反応フロントがまだカラムを通り抜けていません。グラフが水平に安定した時刻 = 反応フロントがカラム全体を通過した時刻です。

pH
 ▲
 │         ___________  ← 入口:早期に安定
 │        /
 │       /
 │      /
 │_____/               ← 出口:遅れて安定(これが「ラグ」)
 └─────────────────▶ 時間 [h]

② 入口と出口の pH 差(ΔpH)

この差が大きいほど、カラム内で活発な反応が起きています。今回のシミュレーションでは、方解石(石灰石)が熱水に溶けることで水がアルカリ性に傾くため、出口の pH は入口より高くなる傾向があります。

ノートなぜ方解石が溶けると pH が上がるの?

方解石(CaCO₃)が溶けると、水中に炭酸水素イオン(HCO₃⁻)が放出されます。このイオンが水をアルカリ性に傾けるため、pH が上昇します。化学的に言うと「緩衝作用(バッファリング)」と呼ばれる現象です。


グラフ 2:Ca²⁺ 濃度と飽和指数(SI)

飽和指数(SI)とは?

飽和指数(Saturation Index)は、ある鉱物が「今の水に溶けやすい状態か、沈殿しやすい状態か」を数値1つで表した指標です。

\[SI = \log_{10} \frac{IAP}{K_{sp}}\]

  • IAP(イオン活量積):今の水にどれだけのイオンが溶けているかを表す数値
  • \(K_{sp}\)(溶解度積):その温度・圧力で鉱物が溶ける「限界量」
SI の値 意味
SI < 0 不飽和 → 鉱物がまだ溶ける余地がある。溶解が進む
SI = 0 平衡状態 → 溶解も沈殿もない
SI > 0 過飽和 → 溶けすぎている。沈殿が起きようとしている

2段グラフの読み方

上段に Ca²⁺ 濃度、下段に飽和指数を縦に並べて描きます。2つを重ねて見ることで、「Ca²⁺ が増えたとき SI はどう変化するか」という関係が読み取れます。

期待される挙動:方解石の溶解が進む場合

  • Ca²⁺ 濃度は時間とともに増加する(方解石から Ca²⁺ が放出されるため)
  • SI はマイナス → ゼロに向かって上昇する(溶解が進むにつれ平衡に近づくため)
  • 入口より出口のほうが、変化が遅れて現れる
ノートSI のグラフで「溶解ゾーン」と「沈殿ゾーン」を色分けする

SI = 0 のラインを境に、ゼロ以下の領域(不飽和=溶解)を青く、ゼロ以上の領域(過飽和=沈殿)をアンバーで塗り分けると、状態の変化が視覚的に一目でわかります。


グラフ 3:鉱物量の空間プロファイル

「時間方向」から「空間方向」へ視点を変える

破過曲線は「時間の流れ」を見るグラフでした。空間プロファイルは視点を変えて、「ある時刻でカラム全体をX線で透視したらどう見えるか」を表します。

横軸にカラムの位置(入口=左、出口=右)、縦軸に各セルの方解石残存量をとり、複数の時刻のスナップショットを重ねて描きます。

3つのパターン

プロファイルの形は流速・反応速度・カラム長のバランスによって、大きく3つに分類できます。

パターン A:前進型溶解フロント
  mol  ▲
       │████████████╲____
       │              ↑ここが溶解フロント
       └──────────────────▶ 位置(入口→出口)

  → 入口側から順番に溶けていく。流速が速く反応速度も速い場合に現れる。
パターン B:均一溶解
  mol  ▲
       │╲_______________
       │ 全体が同じペースで減っていく
       └──────────────────▶ 位置(入口→出口)

  → カラム全体で同時に反応が進む。流速が遅く、流体が十分に拡散する場合に現れる。
パターン C:出口側優先溶解
  mol  ▲
       │____╲███████████
       │         ↑出口側が先に溶ける
       └──────────────────▶ 位置(入口→出口)

  → 入口で反応した流体が出口側に到達してから溶解する。反応速度が遅い場合に現れる。

時間が経つにつれて線が全体的に下がっていれば、シミュレーション全体を通して溶解が着実に進んでいる証拠です。


グラフ 4:4パネルサマリー

全体像を一枚で俯瞰する

個別のグラフを並べて一枚にまとめた「サマリー図」です。論文・レポートの図として使いやすく、全体の傾向を一目で確認できます。

パネル位置 内容
左上 pH 破過曲線(入口 vs 出口)
右上 Ca²⁺ 濃度の時間変化(入口 vs 出口)
左下 飽和指数(SI)の時間変化(入口 vs 出口)
右下 最終時刻における方解石量の空間分布(棒グラフ)

右下の棒グラフが特に重要で、「24時間後にカラムのどこで方解石がどれだけ残っているか」を一目で示します。初期値(0.5 mol)を基準線として引いておくと、どのセルで最も溶解が進んだかが視覚的にわかります。


結果の読み方:3つの着眼点

① 定常状態に達しているか?

破過曲線の出口 pH がまだ変化し続けているなら、シミュレーション時間が短すぎます。反応フロントがカラムを通り切るまで時間を延ばす必要があります。逆に、早々に平衡に達しているなら、より長いカラムや速い流速で実験条件を変えてみると面白い結果が得られます。

② 入口と出口の差は物理的に妥当か?

方解石の溶解が正しく計算されているなら、次の2つが同時に成立するはずです。

  • 入口 pH < 出口 pH(出口のほうがアルカリ性に傾いている)
  • 入口 SI < 出口 SI(出口のほうが平衡に近い)

これが逆になっている場合は、反応速度の設定(速度定数や初期鉱物量)に問題がある可能性があります。

③ 数値的なアーティファクトはないか?

「アーティファクト」とは、実際の物理現象ではなく計算の誤差が原因で現れる異常な値のことです。具体的には以下のような症状として現れます。

  • pH が短時間で激しく上下に振動する
  • Ca²⁺ 濃度が突然マイナスになる

これらは数値不安定のサインです。時間刻みを小さくするか、反応計算のサブステップ数を増やすことで改善できます。


まとめ:4つのグラフで何がわかるか

グラフ 主な着目点
pH 破過曲線 反応フロントの通過時刻、入口と出口のラグ
Ca²⁺ と SI 溶解の速度と平衡への到達、溶解・沈殿ゾーンの判定
鉱物量の空間プロファイル 溶解フロントの形状(A・B・C パターン)、最も溶解が進んだ場所
4パネルサマリー 全指標の時空間変化を一枚で俯瞰

シミュレーションは「動かすこと」よりも「読むこと」に本当の価値があります。4つのグラフを組み合わせることで、カラムの中で起きていた24時間の化学的なドラマが、はじめて「見える」ようになります。


このシリーズの他の記事: